{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<KeysViewHDF5 ['all_st_wp']>\n",
      "<KeysViewHDF5 ['all_st_ep']>\n",
      "<KeysViewHDF5 ['all_st_al']>\n",
      "<KeysViewHDF5 ['all_st_io']>\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEECAYAAAACvbKkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXxU1fn48c+9d5bsCSGBEEJCIiCgUkRArQh1xWLdEE1AsdYFpYIFqyI7CBWoCm3hiyJ1DaIsYn8utNVSLQUVLRXREES2sIQkkBCSSTLJzNzz+2OSgRSISSaTScjzfr3mNXfOPffOMxPIk3PPcjWllEIIIYRoJD3YAQghhGjdJJEIIYTwiyQSIYQQfpFEIoQQwi+SSIQQQvjFEuwAmtu2bduw2+2NPr6ystKv4wNN4vOPxOcfic8/LTm+yspK+vbte8Z9bS6R2O12evXq1ejjs7Oz/To+0CQ+/0h8/pH4/NOS48vOzj7rPrm0JYQQwi+SSIQQQvhFEokQQgi/tLk+EiFEcLlcLg4dOoTT6QzKe9d1rT/YWkJ8ISEhJCUlYbVa632MJBIhRLM6dOgQkZGRdO3aFU3TmvW9KyoqCA0Nbdb3bIhgx6eUorCwkEOHDpGamlrv4+TSlhCiWTmdTtq3b9/sSUT8OE3TaN++fYNbi5JIhBDNTpJIy9WYn40kknoyTcXq/xzE5ZFV94UQ4lSSSOppz1EHT67dzleHy4MdihAiyIqLi3n//feDHcYZbdy4kVWrVgGwaNEihg8fzmuvvcaSJUsC9p7S2V5P0aHeEQxF5Z4gRyKECLbvv/+ef/7zn9x0003BDuU0gwcP9m2vX7+ed999l4iIiIC+pySSeooNt6FpUOx0BzsUIc4Z72w9xOr/HGzSc97Zvwu3X5J01v233XYbf/7zn4mKiuLSSy9lxYoV9O7dm9tuu41bb72V9evXo2kaw4YN45577uGjjz5i+fLlWCwWOnfuzO9//3tefPFFdu7cyapVqxg8eDDTp0/3rZM1Z84cPB4PY8eOJSYmhsGDB7Nx40Z69uzJDz/8gMPh4I9//COdO3c+Y3yLFy9m7969FBYWUlJSwrRp0+jfvz8rVqzgo48+wu12ExkZyeLFizFNk8mTJ5Obm4vL5WL69Ons27ePvXv3EhISQl5eHg899BBjxozhL3/5C4sWLWLNmjW89dZbmKbJNddcw/jx4/3+zuXSVj1ZDJ12YTaOV0iLRIjW7JprruHf//43W7duJSkpic2bN7N7926Sk5P529/+xsqVK1m5ciX/+Mc/2Lt3Lx988AH33nsvb731FoMGDcLhcPDwww9z2WWXkZ6ezoIFCxg9ejSZmZncf//9PPfccwAcPXqUl19+mQcffBCAPn368Nprr3HFFVfw4Ycf1hljSEgIb7zxBs8++yxPP/00pmlSXFzMa6+9xsqVK3G73Xz77be8/fbbdO7cmVWrVjF//ny++eYb3znGjRtHfHw8r7zyCiEhIQAUFhayfPlyVq5cybp16ygtLaWsrMzv71RaJA0QF2Gj2CmJRIimcvslSXW2HgLh+uuv58UXX6RTp05MnDiRzMxMlFIMHTqUBQsWcO+99wJw4sQJDhw4wOTJk1m2bBlvvfUWaWlpXHvttbXOt2vXLpYtW8af//xnlFK+iXxJSUnYbDZfvd69ewOQkJDAsWPH6ozxsssuA6B79+4cO3YMXdexWq089thjhIWFkZeXh9vtZu/evb5LWT169KBHjx6sW7furOc9ePAg3bt39yWWKVOmNOCbOztpkTRAXISdYmmRCNGq9ejRg0OHDrF9+3aGDBlCeXk5GzZsIC0tjW7duvHGG2+QmZnJ8OHD6dGjB6tWrWL8+PGsWLECgI8//hhd1zFNE4C0tDQef/xxMjMzmT17NkOHDgVA1xv/6zUrKwvwJqmOHTuyc+dO/vGPf/CHP/yB6dOnY5omSinOO+88vv32W8CbJH7729/Wed7k5GT27t1LVVUVAI8++ij5+fmNjrOGtEgaIC7Czr6CE8EOQwjhpwEDBnDo0CF0XWfAgAHs3r2bnj17cvnllzNy5Eiqqqro06cPHTt2pE+fPvzqV78iJiaG8PBwfvazn1FVVcWuXbt47bXXmDRpErNmzaKyshKn08nUqVP9ji87O5tf/vKXVFRUMGfOHFJSUggNDWX48OHYbDbi4+MpKCggIyODKVOmcPfdd+PxeJgyZQo//PDDWc8bGxvLgw8+yN13342maVx11VV07NjR73g1pVRAJkYsW7aMf/7zn7hcLkaOHMnAgQN56qmn0DSN7t27M3PmTHRdZ8mSJXz66adYLBamTJlCnz59yMnJ8bvu2fiz3v+cD3bw5hf72Tl3WGO/loBryfczAInPX+dCfMH8DMFeguTHVFRU8Oc//5m4uDhGjhwZtDjO9DOq6+cWkBbJli1b+Prrr3nrrbeoqKjglVdeYd68eUyYMIFLL72UGTNmsGHDBhITE/nyyy9Zs2YNR44cYfz48bzzzjt+173uuusC8bGIi7DjdCvKq9yE2aQxJ4RovHHjxnHiRO0rHGFhYVx44YVBiqjxAvLbcNOmTfTo0YNHHnkEh8PBk08+yerVqxk4cCDgHee8efNmUlNTGTRoEJqmkZiYiMfjoaioiKysLL/qBi6ReDvOjpVWkdxeEokQovHONEGwpbeYziYgvw2PHz9Obm4uL774IocOHWLs2LEopXxruISHh1NaWorD4SAmJsZ3XE25v3XrUllZ2ehlmiuOe2e1/+e77ynrENKocwSa0+kM+jLUdZH4/HMuxOdyuaioqGimiGpTSgXtveujpcTX0OXsA5JIYmJiSEtLw2azkZaWht1uJy8vz7e/rKyMqKgoIiIiao1hLisrIzIyslYfR2Pq1sWfe7Z7ok7AhjzC4zrRq1dCo84RaOfCNfRgkvj8U98+kmD91d3S/+JvKfFZrdYz9pGcTUCG/15yySX8+9//RilFfn4+FRUVXH755WzZsgXwrgXTv39/+vXrx6ZNmzBNk9zcXEzTJDY2lt69e/tVN1DiIuwAHHNUBuw9hBCitQlIi+Sqq67iq6++YsSIESilmDFjBklJSUyfPp2FCxeSlpbG0KFDMQyD/v37k56ejmmazJgxA4BJkyb5VTdQ2lf3kRwtlUQihBA+qo3ZsWOHX8dfNGO9mvru9iaKpun5+/kCTeLzz7kQXzA/Q3l5ud/nePvtt1VVVVUTRHM6f+N75JFHlFJKffPNN2rYsGHqueeeUxMmTFCVlZUNOs+ZfkZ1/dxk6FEDxYQYHCutCnYYQoggWbZsGbfeemuwwzijmpFgmzZtIiMjg9GjRzfL+0oiaaB2oYb0kQjRVLa9BV+vaNpzXnw39D37ZD6Xy8XMmTPJycnBNE0mTJjA8ePHefPNN311/vjHPwIwYcIElFK4XC5mz57N9u3bOXr0KBMnTmTp0qU8//zzfPXVVyiluPfee/n5z3/O6NGjadeuHSUlJdx4441s2rQJp9PJgQMHePDBBxk+fPgZ4zp06BDjx4+nY8eO5OfnM3jwYCZOnMiuXbuYP38+pmn6VgPu16/fGVfxveKKK3jhhRdYu3YtVquVhIQE5s2bx1//+leOHDnCtGnTcLlchISEsGjRImJjY5vkK5dE0kAxIQY50kciRKu1Zs0a2rVrxzPPPMPx48e5++67uemmm3jppZcIDQ1lxowZbNq0iaioKCIjI3n++efZvXs3DoeDO+64gxdeeIFFixbxr3/9i0OHDvH2229TWVnJnXfeyRVXXAHATTfdxHXXXce6detwOBy8/PLL7N+/n4cffvisiQTgyJEjvPrqq0RGRjJq1CiysrLIyclh0qRJnH/++bz//vusW7eOlJQUli9fznvvvYfNZmP+/Pm+Ua19+vThtttuIy4ujuuuu4558+YBsGDBAsaMGcPgwYNZv349O3bsYNCgQU3ynUoiaaCYUIOtR+QuiUI0ib4j62w9BMKuXbvYunUr27dvB8DtdqPrOpMmTSI8PJy9e/fSt29fBg8ezP79+/n1r3+NxWJh7Nixp50nKyvLd/nI7XaTm5sLQGpqqq9ez549AejUqZNvscSz6dGjh2++XJ8+fdi3bx8JCQksXbqUkJAQysrKiIiIaNQqvvv27ePiiy8GYNiwpl3mSVb/baDYUAuOSjcVVbIKsBCtUVpaGjfeeCOZmZksX76cIUOG8MYbb7Bo0SLmzp2L3W5HKcWWLVvo0KEDr7zyCmPHjmXhwoUAaJqGaZqkpaVx6aWXkpmZyeuvv87Pf/5zkpKSfHVqnLr9Y/bu3UtFRQUej4ft27fTrVs3fve73/Hoo4+yYMECevTogVKqUav4nrpS8HvvvUdmZmaDvre6SIukgWJCDMA7l6RLbFiQoxFCNFRGRgbTpk3j7rvvxuFwMGrUKC655BJuu+02wsLCiIqKoqCggKuvvpqJEyfy+uuvo+s6jzzyCAD9+/dnzJgxvPHGG3z55ZeMGjWK8vJyrr32Wr9vaWu1WvnNb37DsWPHuOGGG+jZsyc333wzv/71r2nfvj0JCQkcP368Uav4Pvnkk8yYMYMXXniBkJAQnn32Wb9iraVBY8LOAf4OPXzto/+olEkfqK05RU0UUdM6F4aHBpPE55+2MPw3UA4ePKhuv/32YIehlJLhvwHXLrS6RSId7kKIRli1ahUffPDBaeWPPfZYEKJpGpJIGujkpS2ZSyKEaLj09HTS09PPuK8p+y2ak3S2N9CpfSRCCCEkkTSY1dCICbPKeltCCFFNEkkjxEXYpUUihBDVJJE0QlyETRKJEEJUk0RST8rj4dhLy6GysrpFIp3tQggBkkjqrWr/fo4uXAj/+Y83kUgfiRCt3ksvvcSgQYOorPT+f/7+++/56quvALj66qt95fVVWVnJmjVrAFi3bh0bNmxocExffPEF99xzD6NHjyYjI4NXX30VpVSDz9OcZPhvPVni4rwbxwqJP89OaaUbp8tDiNUIbmBCtGLv7XmPd394t0nPeVv327j5vJvrVff9999n2LBhfPjhhwwfPpyPPvqIuLg4BgwY0Kj3Pnr0KGvWrOGOO+6oc3HGs9m9ezcLFixg2bJldOjQAbfbzaxZs3j55Zd54IEHGhVTc5BEUk96VBRaaCiqsJC4U+6UKMukCNE6bdmyheTkZDIyMnjiiSe44oorePfdd7FarVxwwQW+emdbxv3666+nX79+7Nu3j/bt27N48WJefPFFdu/ezZIlS1BKERcXR3p6OnPnzmX79u24XC7Gjx/Ptddee8aY1qxZw0MPPUSHDh0AsFgsPPXUU9x222088MADfPLJJ757jvTu3ZvZs2fzn//8h0WLFmEYBl26dOHpp5+msrKSqVOnUlpayvHjx7njjjsYNWoUo0ePpmfPnvzwww84HA7++Mc/0rlzZ7+/S0kk9aRpGtaOHakqLCQ+8uS92yWRCNF4N593c71bD02tpuWQlpaGzWYjLy/Pt/x6nz59fPV279592jLu/fr14+DBg7z++ut06tSJjIwMvv32Wx5++GF27drFuHHjWLx4MQAbNmzg+PHjrF27lqNHj7JixYqzJpJDhw6RnJxcqywiIoKKigqqqqqYM2cOa9asoX379ixZsoQjR44wffp0Vq5cSfv27fnDH/7Au+++ywUXXMCNN97I9ddfT35+PqNHj2bUqFGAd1XhqVOnsmjRIj788EPGjBnj93cpiaQBLAkJVBUWEhdRk0ikw12I1ujEiRNs3LiRoqIiMjMzcTgcrFix4rRf4gAdOnQ4bRl3gHbt2tGpUyfAu0T82fpT9u3bR9++fQGIj49n4sSJZ42rQ4cOHD58mN69e/vKHA4HNpuNEydOEBUVRfv27QEYN24chYWFFBQUMGHCBACcTidXXHEFQ4YM4fXXX+ejjz4iIiICt9vtO1/NuRMSEjh27Fi9v7O6SGd7A1g7doRaiUQ63IVojd577z1uv/12XnnlFV5++WVWr17N5s2bfUvEn+pMy7jDmZeH13X9tOPT0tJ8y7eXlpZy//33nzWuO++8kxdeeIGjR48C3rs5/u53vyMjI4P27dtTUlJCcXExAHPnzuXw4cO++5VkZmby8MMPc+mll/LKK6/Qt29fnnvuOW644YaAd9ZLi6QBLAkJUFRErCzcKESrtmbNGn7/+9/7XoeGhnL99ddjGAZvvvkm5513nm/fmZZxP5v27dvjcrl49tlnfTeduuaaa/j8888ZOXIkHo/Htxz9mfTq1YuJEycyceJEPB4Pbreb6667jgceeABd15k5cyYPPfQQuq7Tu3dvLrroIqZOncqYMWNQShEeHs7vf/97NE1j1qxZvP/++8TExGAYxo/eVMsfmmrp48qaWHZ2Nr169WrUscffeou82U/T7V//4pIXvua2izsz+5YLmzhC//jz+ZqDxOefcyG+YH6GiooKQkNDg/Le9dFS4jvTz6iun5u0SBrA0jEBAHd+HnGRdo7KpS0hRANVVVWd8fJWamoqkydPDkJE/pNE0gDWBO8dyFx5ecRHRHKsVDrbhRANY7PZzrpcfEVFRTNH0zQClkhuvfVWIiMjAUhKSiI9PZ3f/e53GIbBoEGDGDduHKZpMmvWLL7//ntsNhtz584lJSWFbdu2+VU3UCwJ1S2SvHziIuPIzi0J2HsJIURrEZBEUjMM7tSse8stt7B48WK6dOnCmDFjyMrK4vDhw1RVVbFq1Sq2bdvG/PnzeeGFF5g5c6ZfdU+dTNSUjHbtwGLBlZ9H/Pn92CiXtoQQIjCJZOfOnVRUVHDffffhdrsZP348VVVVvjHagwYN4vPPP+fo0aNceeWVAPTt25fvvvsOh8Phd91AJRIt/zto397bIrnERqlTlkkRQoiAJJKQkBDuv/9+7rjjDvbv38+DDz5IVFSUb394eDgHDx7E4XD4JvcAGIZxWllj6talsrKS7OzsBn8mq+Mw3T68HS1yICX79uEqLQJgy7YddIhoOV1NTqezUZ+vuUh8/jkX4nO5XEHrC1BKteh+iJYSn8vlatC/s4D8BkxNTSUlJQVN00hNTSUyMtI3iQagrKyMqKgonE4nZWVlvnLTNImIiKhV1pi6dbHb7Y0beuhoDx+CNdJAHS3hou5d4fNjxCR0oVeXmIafL0DOheGhwSTx+ae+w3+DNcS1oqKCr776iiNHjpz1vuk19uzZw6xZs8jMzGTixIksWLAAm80W8PhawvBfq9V6xuG/ZxOQRLJ27Vp27drFrFmzyM/Pp6KigrCwMA4cOECXLl3YtGkT48aNIy8vj08++YRhw4axbds2evToQUREBFar1a+6AREWB4YdS5iiIj+f+HArgNxyVwg/FP/lL5x4Z12TnjP69uHE3HrrWfcPHjy4wedctGiRPyGd8wKSSEaMGMHkyZMZOXIkmqbxzDPPoOs6jz/+OB6Ph0GDBvGTn/yEiy66iM2bN5ORkYFSimeeeQaA2bNn+1U3IHQdohKxhlRR7nIR6/E2P2WZFCFal3Xr1vHvf/+b3NxcEhISOHjwIBdddBGzZ8+moKCAxx9/HKUU8fHxvmOuvvpq/vrXv5KTk1PvlYAN48x9pzt27GDOnDkYhoHdbmfOnDkkJiaydOlSPv74Y0zTZOTIkWRkZJCZmckHH3yApmkMGzaMe+65p0GrEZ8thian2pgdO3Y0/uBXb1RHJwxQO87vqYq3bVcpkz5QizfsarrgmoBfn68ZSHz+ORfiC+ZnKC8vV++8846aMGGCGjhwoCotLVVut1v97Gc/UwUFBWr+/Plq1apVSimlPvzwQ3X33XcrpZS66qqrlNPpVB9++KHauXOnUkqp9957T02dOlUppVTPnj1Vbm6uUkqp9PR09fXXX581httuu833HXz88cdq/PjxKisrS6Wnp6vS0lJVXl6u5syZo3bt2qUyMjKU2+1WHo9HjR49Wu3Zs6dJYvgxZ/oZ1fVzazm9xK1BdBfslu8BC3rhUSJDLLICsBCtVHJysm+wTnx8PJWVlfzwww/ccsstAPTr14+33nqr1jH+rgQMUFBQ4Ot/GDBgAM8//zz79u2jT58+GIZBaGgo06ZNY/369eTm5nLvvfcC3hWLDxw40CQxNDVZ/bchopMI0QuBmtntskyKEK3VmVbvTUtL4+uvvwbwrdh7qoasBHw2HTp0YOfOnQB89dVXdO3albS0NHbs2IFpmrhcLn71q1+RlpZGt27deOONN8jMzGT48OH06NGjSWJoatIiaYjoJCwhHjAM71ySiFRZAViIc8hvfvMbJk6cyPr160lKSjptf0NWAj6buXPnMmfOHJRSGIbBM888Q5cuXbjyyit9rY+RI0fSs2dPLr/8ckaOHElVVRV9+vShY8eOTRJDU5PVfxti9wZYMZzdn1xE2GVX8Ls+6WTnlfDP3/6sSWP0x7kwPDSYJD7/yOq//mkp8cnqv4EU3QUAS0w4rrx84n5qkxaJEOI0ubm5TJo06bTyAQMG8OijjwYhosCSRNIQ0Z0BsEZZcOblERdhp8TpptLtwW6RZVKEqC+lVFCv6QdaYmLiWVf4bekac5FKOtsbwhaO2xaNJczElZ9PXIR3lmuhjNwSot5CQkIoLCwM+O1fRcMppSgsLPTd3bG+pEXSQO6wjljtTpTTSTzeBHLMUUliTPCvawrRGiQlJXHo0CHffcmbk8vlwmq1Nvv71ldLiC8kJOSMAw3qIomkgVzhCViMXADiKrzrh8nsdiHqz2q1kpqaGpT3PhcGK7REcmmrgVxhHbFqxwCIcXhXAJb1toQQbZkkkgZyhSVgsXjvjBh2wptIZHa7EKItk0TSQK6wjlhCTO8ijscKiLRbpEUihGjTJJE0kCssAU0HS7vI6nu326WPRAjRpkkiaSB3WAcALDFhuPPziIuwSSIRQrRpkkgayB3SHnQL1kjDO7s9wi59JEKINk0SSUPpBkQlYgnzVK8AbJM+EiFEmyaJpDGiu2C1VaDKy0mweDhR4aLKbQY7KiGECApJJI0RnYTF8A4BTqj0PheWSatECNE2SSJpjOgk36RE3+z2UuknEUK0TZJIGiM6CUuIy7tZJsukCCHaNkkkjRHdBWuoB4CIkuplUiSRCCHaKEkkjRGdhGaAER2B7bj3EpeM3BJCtFX1Wv3X4XCwfPlyjh49ys9+9jPOP/98UlJSAh1byxVVfYOrmBDU0QIiEi1yaUsI0WbVq0UyZcoUunTpwv79+4mLi2Pq1Kk/ekxhYSFDhgxhz5495OTkMHLkSEaNGsXMmTMxTe9Q2SVLljBixAgyMjLYvn07QJPUDYQqTxXzv5yPw+2AkCgIicYSaeDOq5ndLp3tQoi2qV6JpLi4mBEjRmCxWOjXr9+P3tnM5XIxY8YM31225s2bx4QJE1i5ciVKKTZs2EBWVhZffvkla9asYeHChcyePbtJ6gbKkbIjvJn9JluLt3oLortgDXVX3ynRLvduF0K0WfXuI9mzZw8AeXl56Hrdhy1YsICMjAw6dPCuS5WVlcXAgQMBGDx4MJ999hlbt25l0KBBaJpGYmIiHo+HoqIiv+sGSmJ4Irqmc8R5xFsQnYTFVoFZUkInm5JLW0KINqtefSRTp05lypQp7Nmzh0cffZSZM2eete66deuIjY3lyiuv5KWXXgK89wHWNA2A8PBwSktLcTgcxMTE+I6rKfe37o+prKwkOzu7Ph/7NPG2eA6XHSY7O5uOZjiGVgSEEHU8l/zSiEaftyk5nc4WEcfZSHz+kfj8I/EFRr0Syfnnn8+qVavqdcJ33nkHTdP4/PPPyc7OZtKkSRQVFfn2l5WVERUVRUREBGVlZbXKIyMja7V2GlP3x9jt9kbfyrLboW7kFud6jy+8iDLb+0AIPUM13jxm0q3H+ViN4A6Ea+m36pT4/CPx+Ufia7y6Ely9futdffXVXHPNNb7HLbfccta6b775JitWrCAzM5NevXqxYMECBg8ezJYtWwDYuHEj/fv3p1+/fmzatAnTNMnNzcU0TWJjY+ndu7dfdQMpOSqZPGeet4/olLkkcU7vpMRC6XAXQrRB9WqR/O1vfwO8l6i+++473+v6mjRpEtOnT2fhwoWkpaUxdOhQDMOgf//+pKenY5omM2bMaJK6gZQSlUKFWUGhs5C4qM5YwryJJMZRDHTimKOShOiQgMYghBAtTb0Sic1m821fcsklLFy4sF4nz8zM9G2vWLHitP3jx49n/PjxtcpSU1P9rhsoKVHeuTM5JTnERSehG2BEhhJeUgRWmd0uhGib6pVInn/+eV+ndkFBwY+O2jpXpUR6E8mBkgNckvYT0HQs0SGEFB+DeGQIsBCiTapXIklLS/Nt9+zZkyuvvDJgAbVknSI6YWgGOSU5YFggMhFrhIZZeNSbSKSPRAjRBtWZSDZt2gRAfHx8rfJvvvmGQYMGBS6qFsqiW+hg78CB0gPegugkLKHFVOzPJ6yPIettCSHapDoTyYcffnjWfW0xkQB0CunE/pL93hfRSVith/Ec99ApVJdJiUKINqnORDJv3rwzlhcUFAQkmNagk70TO47twFQmenQSFv1vQBRplEsiEUK0SfXqI/nTn/7EypUrcblcOJ1OunbtWmdr5VzWKaQTTo+TgvICEqKTsIZ4+0W6eErZ5GgX5OiEEKL51Wv41caNG9m4cSM33XQT69evp2PHjoGOq8VKCEkAvCO3iO7im0vSqapUOtuFEG1SvRJJTEwMNpuNsrIyUlJSqKioCHRcLVZiSCIAOaU53j6SUO/S9fEVxRwvr8LlCdxS9kII0RLVK5EkJCSwdu1aQkNDef7553E4HIGOq8WKtcViN+zknPAmEt2q0MPstCsvRikoKpNWiRCibalXIpkzZw6XX345Tz75JB06dGDRokWBjqvF0jWdLpFdvC2SkGiwRWKNtp28d7sMARZCtDH1SiS33347n376KQCjR4+mW7dugYypxUuOTPb2kWiady5JuIa92Hvvdhm5JYRoa+qVSF566SWcTie//OUveeqpp9i6dWug42rRUqJSOFh6EI/p8faThFRhFNYkErm0JYRoW+qVSOLi4rj//vtZvHgxlZWVjB07NtBxtWgpUSm4TBd55XneFonVAUWFWEy3tEiEEG1OveaR/A8bh4YAACAASURBVOUvf+Hdd9/FNE1uv/32s05UbCuSo5IByDmRQ+fozlgtDiCGRHeZLNwohGhz6pVIdu7cyYwZMzjvvPNqlR8+fJjOnTsHJLCWzLecfGkOPz1lLsl5OGQpeSFEm1OvS1tPPfXUaUkEYPLkyU0eUGsQHxpPqCW0elJiEtbqRJJslsmlLSFEm+PXjUWUUk0VR6uiaRrJkcne5eSjk7BU33I3saqEY6XS2S6EaFv8SiQ1N7tqi1KiUrzLyUcmoltBt1uId5ZIi0QI0ea0zVsdNoGUqBQOlR7CpWtoUQlYomy0Kz9OUXkVblkmRQjRhsilrUZKjkrGozzkOnK9/SThEFFS5F0mpVwubwkh2o56JZL33nvvjOWXXXZZkwbTmvhGbtX0k4RUElJcCMgyKUKItqVeiWT16tVnLH/kkUeaNJjWJDnSO5fEN3LLUopRXIRuemR2uxCiTanXPJKqqipuvfVWUlNT0XVv7nn++ecDGlhLFxsSS6Q10nvb3eguWEKq0EyT2MpSmZQohGhT6pVIHn/88UDH0epomkZyVPXijckX++aSxFWckJFbQog2pV6JpEePHmzatAm3241SioKCAgYOHHjW+h6Ph2nTprFv3z4Mw2DevHkopXjqqafQNI3u3bszc+ZMdF1nyZIlfPrpp1gsFqZMmUKfPn3Iycnxu25zSI5KZvvR7d4+kpo7JbpkCLAQom2pVyJ59NFH6dq1K7t27cJutxMaGlpn/U8++QSAt99+my1btvgSyYQJE7j00kuZMWMGGzZsIDExkS+//JI1a9Zw5MgRxo8fzzvvvMO8efP8qnvdddf5/83UQ0pUCn/f/3eqIjpirZ6UmOxxSB+JEKJNqfef7k8//TSpqam8+uqrnDhxos661157LXPmzAEgNzeXuLg4srKyfK2YwYMH89lnn7F161YGDRqEpmkkJibi8XgoKiryu25zSYlKwVQmh9xl6OGhaFaDRFepjNoSQrQp9WqRAFRWVlJRUYGmaZSXl//4iS0WJk2axMcff8yf/vQnPvnkE99M+PDwcEpLS3E4HMTExPiOqSlXSvlV98c+R3Z2dn0/9mmcTqfveOXwzqPZlL2ZzqHx6BEQ6yjkUGGJX+/hj1Pja4kkPv9IfP6R+AKjXonkrrvu4rXXXuOKK65gyJAhXHLJJfU6+YIFC3j88ce58847qaw8+Vd6WVkZUVFRREREUFZWVqs8MjKyVh9HY+rWxW6306tXr3rFfybZ2dm+4xMrE2EHmFEm9vg07BH76FDlwOHCr/fwx6nxtUQSn38kPv9IfI1XV4Kr16WtoUOH8tBDD3HHHXewfv36H71n+1/+8heWLVsGQGhoKJqmceGFF7JlyxYANm7cSP/+/enXrx+bNm3CNE1yc3MxTZPY2Fh69+7tV93mEm2PJsYe471/e3QSVruTiNIiisqq8Jhtd9a/EKJtqbNFkp6eftaFGd9+++2zHnf99dczefJk7rrrLtxuN1OmTOG8885j+vTpLFy4kLS0NIYOHYphGPTv35/09HRM02TGjBkATJo0ya+6zck3BDi6BxZrGaEnilCmSVFZFfGR9maNRQghgqHORLJw4cJGnTQsLIw//vGPp5WvWLHitLLx48czfvz4WmWpqal+120uXaO6suXIFuh2NZYwD5rHQ0ylg2OOSkkkQog2oc5EUnP3wyVLlpy2b9y4cYGJqJVJjkzmvT3vURHRQSYlCiHapHp1tsfFxQHe1X537NiBacoy6TVqFm88YNFJDvV+L3EVJ2QIsBCizahXIsnIyKj1+oEHHghIMK1RclT14o1mJefVtEicxdIiEUK0GfVKJPv27fNtHz16lCNHjgQsoNbGt5x8eS5GbByaoZHgPCGz24UQbUa9EknNCCnwzsN48sknAxZQaxNuDScuNI4DJQfQYpKwRBaS6C7lW7m0JYRoI+qVSDIzMwMdR6uWHJl88gZXYQV0cJZwVC5tCSHaiHolkiuvvJKioiLatWtHcXExNpuNuLg4Zs6cyRVXXBHoGFu8lKgUNh7aCNE/xWr7jNji43JpSwjRZtRrZvuAAQN4//332bRpE+vXr+faa69l+fLlZ5wr0hYlRyVT6CzEEdEBS0gVkaXHOVriDHZYQgjRLOqVSPLy8khLSwMgOTmZI0eOkJKSgmEYAQ2utfANAbbbsYZ5MDxu3MeLZJkUIUSbUK9LW/Hx8Tz33HNcfPHFfP3118TFxbF582asVmug42sVfIlEh6Tq+5K0Lz/B8fIq4iJkdrsQ4txWrxbJ73//ezp06MDGjRvp1KkT8+fPJywsrNFLqJxrukR2AWC/WXHK7HaZSyKEaBvq1SKx2Wz07dvXt7zx9u3bGTBgQEADa01CLaF0DOvIAedRLBHey31xzhMcK62ChCAHJ4QQAVavRDJu3DiOHz9Op06dfDeSkkRSW0pUCjmlB7B07ARapay3JYRoM+qVSAoLC+tcNl54E8nHOR+jteuCJWIfcRXFst6WEKJNqFcfSWpqKvn5+YGOpVVLiUqhuLKYE5EJWMM8xDtLpEUihGgT6tUi2bp1K1dddRWxsbG+sk2bNgUsqNYoOdK7eGNOWCSxtgo6njjBl5JIhBBtQL0SyUcffRToOFo93+KNVisdwzy0P3KcY3JpSwjRBtSZSJYuXcqvf/1rHnvssdNuufv8888HNLDWJikyCV3TOaB5+GmYB6vbRVnh8WCHJYQQAVdnH8nVV18NwFVXXUW/fv0YMGAA27Zt46KLLmqW4FoTm2GjU3gncjzlWKsnJaqCgiBHJYQQgVdnIunZsycA69at47zzzuOzzz7jscceY8OGDc0SXGuTEpXCAecxLNWTEi1FRzFlmRQhxDmuXqO23G43AwYMoKSkhBtvvFFutXsWyZHJ5JQexBIbBUBseTHHy2UVYCHEua1eicTlcjFv3jz69+/PF198gcfjCXRcrVJKVAoOl4OShE6gUT0pURKJEOLcVq9EMn/+fFJTUxkzZgxFRUU8++yzgY6rVaq5f/vBmHiMUE1mtwsh2oR6Df/t2rUrXbt2BWDYsGGBjKdV6xrVFYCckHB+EuoizikLNwohzn31SiQN4XK5mDJlCocPH6aqqoqxY8fSrVs3nnrqKTRNo3v37sycORNd11myZAmffvopFouFKVOm0KdPH3JycvyuGyyJEYlYNAsHLDr9Q110OHGcAzKXRAhxjmvy37rvvfceMTExrFy5kuXLlzNnzhzmzZvHhAkTWLlyJUopNmzYQFZWFl9++SVr1qxh4cKFzJ49G8DvusFk0S10juzMflxYwjzEVZyQe7cLIc55Td4iueGGGxg6dKjvtWEYZGVlMXDgQAAGDx7M5s2bSU1NZdCgQWiaRmJiIh6Ph6KiIr/rXnfddXXGV1lZSXZ2dqM/n9PprPP4WC2W3WW5WEM9hLqrOLgnBz/ersnjCzaJzz8Sn38kvsBo8kQSHh4OgMPh4NFHH2XChAksWLDANzM+PDyc0tJSHA4HMTExtY4rLS31LVPf2Lo/xm63++6r0hjZ2dl1Hn9B6QW8sysbS5h3iHRIhdOv92uoH4sv2CQ+/0h8/pH4Gq+uBBeQDoUjR45wzz33cMstt3DTTTfV6rcoKysjKiqKiIgIysrKapVHRkb6XTfYUqJSqPA4cUR6YzMLZNVkIcS5rckTybFjx7jvvvt44oknGDFiBAC9e/dmy5YtAGzcuJH+/fvTr18/Nm3ahGma5ObmYpomsbGxftcNtpohwEcSvCslG4VHgxmOEEIEXJNf2nrxxRcpKSlh6dKlLF26FICpU6cyd+5cFi5cSFpaGkOHDsUwDPr37096ejqmaTJjxgwAJk2axPTp0xtdN9h8qwDHx9Kbo4QcL8Q0Fbqu/ciRQgjROjV5Ipk2bRrTpk07rXzFihWnlY0fP57x48fXKktNTfW7bjAlhCVg1a3khIdyYYgitqKY4goXseG2YIcmhBABEbxJF+coQze8a24ZGtZQNx0qjsukRCHEOU0SSQAkRyVzQDmxhrnpWFEkN7gSQpzTJJEEQEpUCgeqTmAJNWlfUSKTEoUQ5zRJJAGQHJVMlXJTEWES5qqkqLAk2CEJIUTASCIJgJRI78itokjv67LDuUGMRgghAksSSQDUDAHOiw0FoPJIXjDDEUKIgJJEEgAdwjoQagklJy4CkNntQohzmySSANA0jS6RXdgVawfAOCaz24UQ5y5JJAGSEpXCXrtCtynsxwuDHY4QQgSMJJIASY5M5pCnAkuYm+jSoyilgh2SEEIEhCSSAEmJSsGNwhNu0qG8kBMVrmCHJIQQASGJJEBqRm6VRShinSUcldntQohzlCSSAKlZTr4oEsIrKzlaJJMShRDnJkkkAdI+pD3h1nCORBkAFB88EuSIhBAiMCSRBIimaSRHJrO/eghw2SGZ3S6EODdJIgmgrlFd+aGdt0VSJbPbhRDnKEkkAZQclczOSBMAj8xuF0KcoySRBFBKVArldg2sJvpRSSRCiHOTJJIAqhm55Q5XRByXS1tCiHOTJJIA6hrVFQBHhCKqRNbbEkKcmySRBFC0PZpoawTHIyGqrESWSRFCnJMkkQRYSmQyR6J0wpxOTpRUBDscIYRocpJIAiw5Jo2cGB0NOHpA5pIIIc49AUsk33zzDaNHjwYgJyeHkSNHMmrUKGbOnIlpeofELlmyhBEjRpCRkcH27dubrG5LkhyVzN5o79d8fP+hIEcjhBBNLyCJZPny5UybNo3KSu9ChfPmzWPChAmsXLkSpRQbNmwgKyuLL7/8kjVr1rBw4UJmz57dJHVbmq5RXSmK0gC5d7sQ4twUkESSnJzM4sWLfa+zsrIYOHAgAIMHD+azzz5j69atDBo0CE3TSExMxOPxUFRU5HfdliY5KpnCSO+2M1fW2xJCnHssgTjp0KFDOXTo5GUcpRSa5v2rPDw8nNLSUhwOBzExMb46NeX+1v0xlZWVZGdnN/qzOZ3OBh3v9Dgpt4PHoqjYv8ev967X+zUwvuYm8flH4vOPxBcYAUkk/0vXTzZ8ysrKiIqKIiIigrKyslrlkZGRftf9MXa7nV69ejX6s2RnZzf4+NhvIiiLOE50aaFf710fjYmvOUl8/pH4/CPxNV5dCa5ZRm317t2bLVu2ALBx40b69+9Pv3792LRpE6Zpkpubi2maxMbG+l23JUoJS6AoEkKLZVKiEOLc0ywtkkmTJjF9+nQWLlxIWloaQ4cOxTAM+vfvT3p6OqZpMmPGjCap2xKlRKeSF/U9CfuLgx2KEEI0uYAlkqSkJFavXg1AamoqK1asOK3O+PHjGT9+fK2ypqjb0qS078nh6L9zaXkFyu1GszRL/hZCiGYhExKbQXJ0KkWRGpoC17FjwQ5HCCGalCSSZpASleIbAnziwOHgBiOEEE1MEkkz6BLZhdz23mHKhUtfwHQ6gxyREEI0HUkkzSDMGoYrNoTPf+bG3PIZBx8cg8fhCHZYQgjRJCSRNJPkkPZ8crHGGz8dQdnWrRz41X24jx8PdlhCCOE3SSTNJC0yiRyrhfJeycwe8Esc2TvZd/doXAUFwQ5NCCH8IomkmaTE9eK4YfB0zFtcdfvlTL/sPhw5B/k+fRRVh6QDXgjRekkiaSYpCf0A2F+8m4d33suzD3Rj8Y0TKD9WxLfD7+TE97uCHKEQQjSOJJJmckH7C7Abdh5LSuGLqHZ0//tolv38GF+Mm4Oz0sX36Xfx3398EewwhRCiwSSRNJOO4R15c9ibRITEMCa0ij90HwCfzuXRyDewL3yeKsOKZ+JYlv3fOpwuT7DDFUKIepNE0ozOjz2ft298m+Hdh/Oy6wi/7DWAg/s/pf+3j9DnxZl4otsxcOksnpj0Et8eOhHscIUQol4kkTSzMGsYs346i+eHPM9+s5w7kpNZb/EQ9bfR9J9zB0ZyCg/+bSkLpr7Ioo934fK0vNsHCyHEqSSRBMn1Xa9n7U1r6dG+F5MidaYld6Nq42QuGN2JsN49mfzVG3z36lvctnQz3+f9+A27hBAiWCSRBFFiRCKvDH2Fh3/yMO9rZdzZrTff719H2pW5RPS7iMf/+zYXbPmImxZv4sV/7cFjqmCHLIQQp5H1zIPMolt4pO8jDEwYyOR/T+aupCQmnCjgrp4HOaIP5N4v19I9QmfuX00+ysojfUAX+iW347z4CHRdC3b4QgghiaSlGJAwgLU3rWXmZzN5Tv2TzyMimePegG4O5IpPV/PmzRq/KbyCSe98C0BkiIW+XWLol9yOfint6NslhuhQa5A/hRCiLZJE0oLEhMTwh6v+wOrvV/Psf57ljrQ0ngnZRqqnK7y3ivfuVJTfPoLtDo2vCj18fegEi//5AzVXvLp1iCAtSuPq0nAuTm5H9w7SahFCBJ4kkhZG0zTSe6bTr2M/ntz4JA+pKn4VdoyRHivFq1fD6tX0AnrpOr9q1w4tph0VYZEUWcM4fCCUPZUG//4ojA/sEVRFRNO5ayfO65bEBb26cHHX9sSE2YL9EYUQ5xhJJC1U93bdWXnjSp776jle3bWaL4dqPJNYQPuyUDxVFtxVFjyVx/E4T2AthPAKRacKRT/nmYcLKyDHbuF7m43ykHDKwqJxRMRRERmHio5Bj4nBaNcOW2wsIfHtCY9vT1RsNNFhNqJDrUSHWgmx6miatHCEELVJImnBQi2hTL98Oj9N/CkzPptORj+dXpYoOmo2OmKlIwYd0enogY4mxJkenMXF2E0Nt6MKj6MKd5kLd7kbl8ODp8KJcipvAirOxZ2v46nUQZ05Obg1nWO2cPbYwymxhVNmD8dts2Na7SibDWx2NLsdLcT7bISEYISEYAkNwRoWijXUji0sFFt4GCFhIdgjwjiSl0d5pRXdYsGwGOgWC7phoFutGFYDXdcxNA1D19B1DV0DQ/NuG5qGzaITYjUw5JKdEC2GJJJW4JqUa7gg7gJe/OZFckpy+K48nw1l+VSZVbXqGbpBTMcYOkd3pmNYRzqGdSQhPMG7Hd6RDqHxRJomYY5jWEuPwIlDqOM5mAX78eQfwl1wBM+JE3gqdd+jfaWDKncoVVU23OU6qkSB2wSPQnObaA0ckhx9lnIFuKufTU3H1DVMTcPUdFT1s6lpuAwLlbqVSouNSksIVRY7HqsN02rzPttOJjnsIeBLdN4kp9ts6HYbhs2GYbN6n+02LNWPgqJCjpXoWO02byIMsWOxWbBbDKyG7ktuGt5nNNA1DY3qZ817eVI7YzlYdV36rcQ5RxJJK5EQnsCsn87yvVZKUVxZTEF5Afnl+eSV5ZFfns+u3F04LU52F+9m0+FNVLgrzng+i24hzBJGmDWM0KhQwmITCevTjVDdRpjSCDPdhLqrCKtyElpVRpizhNCKYiweN7rpwaJMdNOD4a5+VHkwPArDDRa3hu7W0D0Kw62heUD3eMsspsJigm6CxQRDKXRTAwVKgfJtnyzD1Kr3efcrj4bpqX6u0vFUGHg8enU54AHcisb8uk46Q5mJRqlu4NYNPJqOqemnPGuYmlH9rOPRdUyqn2vKTjnGo+ko3cDUdZRhQRmG96EbYBgow4Jm8W5jWNAsFrBY0AwDzWLB6XLxeUQEmmGAxUA3LGAx0KqfdYsFrfocusWKbvEep1u9dXTDQNM1NF33PXRdh+oEp+k6hlH92jC8CdF3jIHFYniTrtWC1WbFYrdhtejYDB2LoVFc4eFEuQuLoWE1dKyGJpdD2wBJJK2Upmm0C2lHu5B2nB97vq88OySbXr16Ad5k43A5yC/LJ788n4LyAhwuBxXuCspd5ZS7y0/bPuos8pa5y6lweZ89ygMaEAZgVD+all23YtOt2HTLyW2tZtvifWgGhlIYpgfNdKOb3qSme1zebY/r5MNdhVFVha2qCotLYXFrWN1g9YDNA1aPwuLxvrZ6wFLrodA9mjfheTR0U0czdeymjqY0b+JDQzM1NKWhKUB5X3PK69qJETCrk6VbgQmYCmXiLax+7X3ULmvJv4YV4ESjTPcmSaums7M64XpqnjUDT3XyNDUdNM2bSDUNhbf1qTS9+qH5WqBUbyvdQOE9hppWqu6tX3MeUzNQNefRNdQpr03dqD6Pjsvt4R82m7d56Gs9es+raYCmo+k1+/TqMu3kdnXS1fTq46qTMYZe67V2aj1dh5pkrFXPAdeqE2z1ZwJA1yguPsG2Tbt8sdU0ZbXq70LXdAxDQzMMb+u4+g8BQ9fRjerXRvVrXcM4pZ5h0bnwkp6E2Jp+moAkknOYpmlE2iKJtEXSrV23Rp1DKYXLdPmSjcf04FEeTGXiVm5MZeJRHjymp9Z2TZ1T97mVm5yDOcQlxFHpqaTKU0WVp8q7bZ6y7fmfbbOKMk8lxz1VePCgNIWpm5ga3ktghgWFgUdZMJUdpRQe5UEphalMzOpYTOXBbXpwKXcTf9O1vjEADBSGAgsaFvA9NOXNyTr4EoSO9xKYhla9X/nqAOjKm/QMD+imiaWmNefxvofuUSdfmzUPMKpbfsYpD115z6crbyyG0mq91gHdVN6ECdXlqvq4s5zzlOf/3baY1ec1qU68+JKvN+HiK0OBZp4sp6Ycbzm+cs13HGd8aNXXSE85T+B+4K3KJ1ffwM+XLmry87b6RGKaJrNmzeL777/HZrMxd+5cUlJSgh3WOUPTNGyGDZthI4YYv8+XXZFNr+69miCyxlNK4VZuXB4XLtP7cJtuXKaL73/4ni5du/he19rvceFSLjymB7fpxm268SiPb79HnSx3m27c6uS2x/T4Eq9SCoU6+Yw34aHwbavqhPS/9UscJYRHRPjKPShcysR5yvm8+7yP07eVb9ubbE/WMZWJielLwJ6acrzHeKq3PdXbHkzcqpUsKnpKsqxJ5jUJRqvu5vvRck4ef+q59DO8rtk2zOpkr8A4JVn76kGtJK4phY43uZ8a68n3UL7Y9FPirP1a+d6j9jkU/S93BOTrbfWJ5B//+AdVVVWsWrWKbdu2MX/+fF544YVghyVaME3TsGpWrPrpTXxHqKPWpcKWJjv75KXLluLUlueOnTs4r/t5vjK3Wd1qPSWRmqcksrO9rmlR/m+dmiRbK2melkRPJt6axFiznZeXR0JCAuD9g+JUNcn7bE77A6D6nGZ1Mq15bxPvHwWnJuVafyBU/8Fwanw1719YWEhMu5jaf0j8z3vW1K0pqxUbClMpPMrbFFOnHA8Q33Nk0/zQ/0erTyRbt27lyiuvBKBv37589913QY5IiLZF13R0TceqWwk1Qom2n21sXvBlq2x69WxZifhULfEPhfpo9YnE4XAQERHhe20YBm63G4vlzB+tsrKS7OzsRr+f0+n06/hAk/j8I/H5R+LzT0uP72xafSKJiIigrKzM99o0zbMmEQC73e5Xxm/pfzFIfP6R+Pwj8fmnJcdXV4Jr9fcj6devHxs3bgRg27Zt9OjRI8gRCSFE29LqWyTXXXcdmzdvJiMjA6UUzzzzTLBDEkKINqXVJxJd13n66aeDHYYQQrRZrf7SlhBCiOCSRCKEEMIvkkiEEEL4RVP/O73zHLdt2zbsdnuwwxBCiFalsrKSvn37nnFfm0skQgghmpZc2hJCCOEXSSRCCCH8IolECCGEXySRCCGE8IskEiGEEH6RRCKEEMIvrX6trUD4sdv3rl69mrfffhuLxcLYsWO56qqrmjU+l8vFlClTOHz4MFVVVYwdO5ZrrrnGt//VV19l7dq1xMbGAjB79mzS0tKaNcZbb72VyMhIAJKSkpg3b55vX7C/v3Xr1vHuu+8CJ+9Ps3nzZqKiogCYO3cu//3vfwkPDwdg6dKlvs8SaN988w3PPfccmZmZ5OTk8NRTT6FpGt27d2fmzJno+sm//ZxOJ0888QSFhYWEh4ezYMEC38+8OeLLzs5mzpw5GIaBzWZjwYIFxMXF1apf17+DQMeXlZXFww8/TNeuXQEYOXIkw4YN89UN9vc3ceJEjh07BsDhw4f5yU9+wqJFJ++nrpRi8ODBvvj79u3Lb3/724DG12hKnObvf/+7mjRpklJKqa+//lo9/PDDvn0FBQXqF7/4haqsrFQlJSW+7ea0du1aNXfuXKWUUkVFRWrIkCG19v/2t79V3377bbPGdCqn06luueWWM+5rCd/fqWbNmqXefvvtWmUZGRmqsLCw2WN56aWX1C9+8Qt1xx13KKWUeuihh9QXX3yhlFJq+vTp6qOPPqpV/5VXXlF/+tOflFJKffDBB2rOnDnNGt9dd92lduzYoZRS6q233lLPPPNMrfp1/TtojvhWr16tXn755bPWD/b3V6O4uFjdfPPNKj8/v1b5/v371UMPPRTQmJqKXNo6g7pu37t9+3YuvvhibDYbkZGRJCcns3PnzmaN74YbbuA3v/mN77VhGLX2Z2Vl8dJLLzFy5EiWLVvWrLEB7Ny5k4qKCu677z7uuecetm3b5tvXEr6/Gt9++y27d+8mPT3dV2aaJjk5OcyYMYOMjAzWrl3bbPEkJyezePFi3+usrCwGDhwIwODBg/nss89q1T/13+ngwYP5/PPPmzW+hQsX+m7C5PF4Tlsxoq5/B80R33fffcenn37KXXfdxZQpU3A4HLXqB/v7q7F48WLuvvtuOnToUKs8KyuL/Px8Ro8ezYMPPsjevXsDGp8/JJGcwdlu31uz79TLHOHh4af9Aw208PBwIiIicDgcPProo0yYMKHW/htvvJFZs2bx+uuvs3XrVj755JNmjS8kJIT777+fl19+mdmzZ/P444+3qO+vxrJly3jkkUdqlZWXl3P33Xfz7LPP8uc//5mVK1c2W6IbOnRorbt7KqXQNA3wfk+lpaW16p/6XZ5pf6Djq/nF99///pcVK1Zw77331qpf17+D5oivT58+PPnkk7z55pt06dKF//u//6tVP9jfH0BhYSGff/45w4cPP61+fHw8Y8aMITMzk4ceeognnngioPH5QxLJGdR1+97/3VdWVtZsieMWsAAAA/BJREFU189PdeTIEe655x5uueUWbrrpJl+5Uopf/vKXxMbGYrPZGDJkCDt27GjW2FJTU7n55pvRNI3U1FRiYmI4evQo0HK+v5KSEvbu3ctll11Wqzw0NJR77rmH0NBQIiIiuOyyy4LWYjq1P6SsrMzXh1Pj1O/yTPubw/r165k5cyYvvfTSaf0Ldf07aA7XXXcdF154oW/7f/8ftITv729/+xu/+MUvTruqAHDhhRf6+j779+9Pfn4+qoWuaCWJ5Azqun1vnz592Lp1K5WVlZSWlrJnz55mv73vsWPHuO+++3jiiScYMWJErX0Oh4Nf/OIXlJWVoZRiy5Ytvv9MzWXt2rXMnz8fgPz8fBwOB/Hx8UDL+P4AvvrqK37605+eVr5//35GjRqFx+PB5XLx3//+lwsuuKDZ4wPo3bs3W7ZsAWDjxo3079+/1v5+/frxr3/9y7f/kksuadb4/t//+3+sWLGCzMxMunTpctr+uv4dNIf777+f7du3A/D555+f9nMM9vdXE9fgwYPPuG/JkiW8/vrrgPcyYWJioq+F2tLIqK0zONPte1999VWSk5O55pprGD16NKNGjUIpxcSJE5t9NeEXX3yRkpISli5dytKlSwG44447qKioID09nYkTJ3LPPfdgs9m4/PLLGTJkSLPGN2LECCZPnszIkSP5/+3cvWoiYRiG4WcIcwYDErCS8RA8gTSpYqwCthaC4l+joigIYprRVCmmEkKKFAFrGwvxQOYANGLpD0qKJWF3ky3kI44L99UNNi+jcDMf42tZlu7v7/X8/Hw290+SgiBQNBr9vP79+725udHd3Z1s29bt7a3i8fjJ55Oker2udruth4cHxWIxXV9fS5IymYx831c6nVa9Xlc6nZZt2xoMBiebbb/fq9fr6fLyUsViUZKUSCRUKpVUq9VUqVS+/R38fbTzkzqdjrrdrmzbluM46na7ks7j/n0IguBLhD/my2azqlarmk6nuri4+PE33kyw/RcAYISjLQCAEUICADBCSAAARggJAMAIIQEAGCEkQAhGo5H6/f7n9Waz0evrq6RfKzNeXl7CGg04GiEBzsB8Pv8MCfC/4Q+JQIiWy6Xy+byWy6UWi4UeHx8lSZPJROPxWKvVSuVyWVdXVyFPCvwbTyRASN7e3pTL5dRoNDQcDuW6rgqFgiQpEono6elJzWaTYy6cPUIChGQ2m2m73epwOHz57GMvlOM4Wq/Xpx4NOAohAUKSSqXkeZ5arZYsy/ojKOe6nA/4DiEBQuS6rpLJpHzf1263k+d5YY8EHI2ljQAAIzyRAACMEBIAgBFCAgAwQkgAAEYICQDACCEBABghJAAAI+8m3O0iuP0x+QAAAABJRU5ErkJggg==\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# %%\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import h5py\n",
    "import scipy.io\n",
    "from functools import reduce\n",
    "import operator\n",
    "from scipy.sparse.linalg import svds\n",
    "from sklearn.decomposition import TruncatedSVD\n",
    "import seaborn as sns\n",
    "import os\n",
    "import os\n",
    "\n",
    "out_num = 20\n",
    "for a in range(0, 4):\n",
    "    area_name_list = ['western_pacific', 'eastern_pacific', 'Atlantic_Ocean', 'indian_ocean']\n",
    "    file_name_list = ['split_wp', 'split_ep', 'split_al', 'split_io']\n",
    "    variable_name_list = ['all_st_' + _ for _ in ['wp', 'ep', 'al', 'io']]\n",
    "    with h5py.File(file_name_list[a] + '.mat', 'r') as f:\n",
    "        print(f.keys())\n",
    "        raw_variables = f[variable_name_list[a]][()]\n",
    "    variables_reshaped = np.reshape(raw_variables, (reduce(operator.mul, (raw_variables.shape[0:-1])), raw_variables.shape[-1]), 'F')\n",
    "    svd = TruncatedSVD(n_components=out_num)\n",
    "    singular= svd.fit(variables_reshaped).singular_values_\n",
    "    plt.plot(singular)\n",
    "    plt.legend(area_name_list)\n",
    "    plt.xlabel('kth')\n",
    "    plt.ylabel('singular_value')\n",
    "    f.close()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  1   2   3 200]\n",
      " [  5   6 100   8]\n",
      " [  9  10  11 300]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# 创建一个示例矩阵\n",
    "matrix = np.array([[1, 2, 3, 4],\n",
    "                   [5, 6, 7, 8],\n",
    "                   [9, 10, 11, 12]])\n",
    "\n",
    "# 定义要修改的下标列表\n",
    "indices = np.array([[1,0,2],\n",
    "                    [2,3,3]])\n",
    "\n",
    "# 将下标列表的行和列分别存储在两个数组中\n",
    "rows = indices[0, :]\n",
    "cols = indices[1, :] - 1  # 减1，因为Python中的索引从0开始\n",
    "\n",
    "# 定义一个与下标列表数目相同的值列表\n",
    "values = [100, 200,300]\n",
    "\n",
    "# 将值列表中的值赋给对应的下标\n",
    "matrix[rows, indices[1, :]] = values\n",
    "\n",
    "print(matrix)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import h5py\n",
    "import scipy.io\n",
    "from functools import reduce\n",
    "import operator\n",
    "from scipy.sparse.linalg import svds\n",
    "from sklearn.decomposition import TruncatedSVD\n",
    "import seaborn as sns\n",
    "import os\n",
    "import os\n",
    "area_name_list = ['western_pacific', 'eastern_pacific', 'Atlantic_Ocean', 'indian_ocean']\n",
    "singular_list=np.load('singular.npy')\n",
    "for i in range(4):\n",
    "    plt.plot(singular_list[i])\n",
    "plt.legend(area_name_list)\n",
    "plt.xlabel('kth')\n",
    "plt.ylabel('singular_value')\n",
    "plt.show()\n",
    "print(np.sum(np.square(singular_list),axis=1))"
   ],
   "metadata": {
    "collapsed": false
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
